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Abstract We investigate simulated turbulent flow 
within thermally driven stellar convection zones. Dif- 
ferent driving sources are studied, including cooling at 
the top of the convectively unstable region, as occurs 
in surface convection zones; and heating at the base by 
nuclear burning. The transport of enthalpy and kinetic 
energy and the distribution of turbulent kinetic energy 
dissipation are studied. We emphasize the importance 
of global constraints on shaping the quasi-steady flow 
characteristics, and present an analysis of turbulent 
convection which is posed as a boundary value problem 
that can be easily incorporated into standard stellar 
evolution codes for deep, efficient convection. Direct 
comparison is made between the theoretical analysis 
and the simulated flow and very good agreement is 
found. Some common assumptions traditionally used 
to treat quasi-steady turbulent flow in stellar models 
are briefly discussed. The importance and proper treat- 
ment of convective boundaries are indicated. 
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1 Introduction 

While the equations governing the dynamics of non- 
magnctized stellar plasma are well known, a fundamen- 
tal understanding of fully developed turbulent flow re- 
mains elusive. A Reynolds decomposition, whereby the 
properties of the stellar plasma are separated into mean 
and fluctuating components <p = {(f>) + </>' provides some 
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insight into the problem. Decomposing the kinetic en- 
ergy equation (formulated by the product of the veloc- 
ity and the momentum equation) in this way and taking 
temporal and angu lar averages (indicated by the oper- 
ator 0) results in (jMeakin fc Arnettl[2007bl ) 



d t ( P E K ) + V ■ (p E K u ) = 

-V • (F p + F K ) + (p'V ■ u'} + (W 6 > - (e K ) 



(1) 



which is the full non-linear governing equation of inter- 
est. The primary goal of any stellar turbulence theory is 
to model the terms of this equation, including the rate 
of buoyancy work Wj, = p'w! ■ g, the kinetic energy flux 
Fk = u'Ek, the pressure correlation flux F p = u'p', 
the work done by pressure fluctuations p'V • u', and 
the rate at which kinetic energy Ek , is dissipated ck- 
Differential rotation and circulation currents introduce 
additional sources of turbulence and transport terms. 

It is standard practice to ignore or grossly approx- 
imate most of these terms in stellar evolution calcu- 
lations. For instance, mixing length theory (MLT) ig- 
nores ek , F p , Fk, and p'V-u' and approximates the in- 
tegral of Wfe over a mixing length as the product of local 
properties of the flow. The time dependence expressed 
by the left hand side of Eq. Q]is also dropped. Though 
still not widely used, some strides have been made to 
compensate for these deficiencies through embellished 
MLT type algorithms, most notably to address the is- 
sue of time d ependence and the non-local nature of tur- 
bulence (e.g. |Goughlll973 IUnnolfl98lt lEggletonl Il983l 
KuhfusJl98a iDeng et al.ll2006l and references therein) . 



In order to develop a more realistic physical descrip- 
tion of stellar turbulence, which is o f central impor- 



tance to modeling stellar p ulsation fe.g. lBelkacem et al 



2006t ISamadi et ail I2OO 



r p 
)9) 



a better understanding of 
these non-linear terms is needed. A powerful method 
for gaining insight into this physics is analyzing fully 
non-linear simulation data. In the following we present 
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a few select models from a new suite of turbulent stellar 
convection simulations designed to this end and briefly 
discuss the origin of the kinetic energy flux and its rela- 
tionship to the kinetic energy dissipation and the large 
scale topology of the flow. 



2 Simulation Setup 

The initial conditions used in our reactive hydrody- 
namic simulations are based on a 23 M star model 
which has been evolved with th e TYCHO stellar evolu- 
tion code jArnett et al.ll2009bl ) to an age of ~ 2 x fO 6 
yr, at which point oxygen is b urning in a shell that over- 
lies a silicon-sulfur-rich core (jMeakin fc Arnettll2007b , 
2006). Variations in the driving source are made in 



order to study how this impacts the global characteris- 
tics of the flow, and in turn how this affects the trans- 
port terms. Three models are presented including two 
in which heating (by nuclear burning) is present and 
one in which convection is driven by a cooling region 
near the top of the convection zone (similar to radiative 
losses in surface convection zones). The heating and 
cooling profiles are presented in Fig 1. The heated and 
cooled regions (0.44 < r/10 9 cm < 0.64) are initially 
nearly adiabatic, and thus neutrally buoyant, while the 
surrounding layers are stably stratified. The fully com- 
pressible, reactiv e Euler equations are soly ed using the 



PROMPI code (jMeakin fc Arnettl l2007bl ) which is a 



descendant of the PR OMETHEUS piecewise parabolic 
method (PPM) code jFrvxell. Miiller. fc Arnettl 1 19891 ) 



adapted to parallel computing platforms. We work 
within the implicit large eddy simulation (ILES) frame- 



Boris 


2007; 


Asp den 


2007tlBenzi et al. 



20081 ). The sensitivity of our results on resolution are 



tested within limits of computational cost by a series 
of higher resolution runs. A summary of simulation 
properties is presented in Table 1. 



3 Thermal Relaxation 

The time evolution of the kinetic energy is presented 
in Fig. 1. The convective turnover time r c = 2L/v is a 
little less than t c ^100 s for all of the models studied. 
After an initial transient comparable to r c the mod- 
els attain a quasi-steady state. The strong damping 
present in turbulent flow (see e.g. lArnett et al.ll2009al ) 
ensures that this state is reached within ~ r c . Time 
averages for analysis are performed over intervals that 
encompass ~ 2t c and are summarized in Table 1. 

The simulated convection is very efficient in all cases 
and deviates only mildly from an isentropic state. Since 



the simulated regions are not in thermal balance (there 
is either a net heating or cooling) the entropy will 
change over time and the fluxes and flows within the 
convection zones adjust to establish an isentropic state 
at each moment. In this situation the rate of entropy 
change at any one location is equal to that of the mass 
averaged rate over the convection zone, s = (s) , where 
(•) indicates a mass weighted average over the turbu- 
lent region. From the first law (SQ = dU + 5W) and 
the fundamental thermodynamic relationship (dU = 
TdS - SW) 



e„ + € K dLc/dm 



= (*), 



(T) 7 



T T 
and the convective luminosity is 

L c (m) = I (e„ + e K - T(s) m )dm' 

Jm v ' 



(2) 



(3) 



where e n is the local heating or cooling term (see Fig. 1). 

The convective flux found from this relationship is 
compared to the simulation data in Fig 2 for all three 
models. The good agreement shows that thermal relax- 
ation is not necessary to study turbulent convection but 
can be incorporated into the analysis. A much more 
important effect than this slow thermal relaxation is 
the luminosity associated with boundary layer mixing 
events which is as large as ^40% of peak in model cl. 

The kinetic energy dissipation ek is required to ap- 
propriately calculate (s) m and is included in standard 
stellar evolution (i.e., MLT) only implicitly through the 
structure variable V = dlnT/dlnP. The distribution 
of ck throughout a convection zone, however, is in- 
timately related to the resulting kinetic energy flux, 
which we discuss next. 



4 Kinetic Energy Flux and ek 

In quasi-steady states where F p and p'V • u' are not 
important the kinetic energy flux (or luminosity) can 
be found by integrating Eq. 1 



r-Mo+m dr , 

L K (m)= / [Lc^ad-fr - e K dm') 

J M \ Hp I 



(4) 



with dr' = dm! / Airr 12 p. In this expression the ra- 
dial component of the rate of buoyancy work is writ- 
ten in terms of the convective energy flux with Wb = 
Fc^j 'ad/ 'H p for pressure scale height H p and adiabatic 
gradient V a d = (d\nT/d\nP) s . This expression for Wb 
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Table 1 Selected Simulation Model Parameters 



Model ID 


AO, A(f> a 


Zoning 


t 6 


Comments 




[deg.] 


[n r x ng x rty] 


["] 




hi 


30,30 


200x50x50 


[300, 500] 


narrow, static heating profile 


hl.z2 


30,30 


400x100x100 




model hi with moderate resolution increase 


hl.zl 


30,30 


800x200x200 




model hi with high resolution 


h3 


30,30 


200x50x50 


[375, 575] 


broad, static heating profile 


cl 


30,30 


200x50x50 


[200, 400] 


static top cooling profile 



"The computational domain is centered on the equator so that the domain extends A8/2 degrees above and below the equator. 
^Provided is the time interval over which averages are performed. 
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Fig. 1 (left) Heating / cooling profiles for the models listed in Table 1. (right) The time evolution of the total kinetic 
energy in the simulation domain. The two additional high resolution models hl.zl and hl.z2 are indicated by the magenta 
crosses and the orange line, respectively. 
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Fig. 2 Convective flux: (left) time averaged simulation data and (right) calculated from the background structure as 
described by Eq. 3. 
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Fig. 3 Kinetic energy profiles: (left) total value of Ek and (right) effective isotropic value Ek,iso described in §4. The 
line segments indicate the two assumed distributions also described in §4. 




Fig. 4 Kinetic energy flux: (left) time averaged simulation data and (right) values calculated directly from background 
structure using Eq. 4. The line thickness in the theoretical profiles indicate the assumed distribution of Ek,iso shown in 
Fig. 3: the (thin line) is the uniform case and the (thick line) the case with a gradient. 
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can be calculated directly from the background struc- 
ture (Eq. 3) and is appropriate for small density fluctu- 
ations that can be linearly related to temperature fluc- 
tuations using the isobaric thermodynamic derivative, 
a good approximation in most cases of deep, nearly adi- 
abatic convection. From Eq. we see that the kinetic 
energy flux is the residual between buoyancy driving and 
viscous dissipation. 

Globally, the integrated dissipation J eKdm is con- 
strained by both the thermal state evolution, Ts (Eq. 2, 
3), and the balance with buoyancy driving (Eq. 4, not- 
ing that Lxirtop) = £#r(n>ot) = 0). The radial profile 
of tK is determined by th e topology of the convective 
flow. lArnett et al.l (|2009al ) show that the dissipation is 
well described by the properties of the isotropic compo- 



nent of turbulence, vf 



2 

iso 



m + vi) with e K ~vtJl d 



where Id is the largest scale of motion in the flow and 
ve and Vtf, are the non-radial velocity fluctuations. 

In Fig. 3 we present the radial distribution of the ki- 
netic energy from the simulation data. The first panel 
shows the total Ek and the second panel shows the 
horizontal component scaled to an equivalent isotropic 
value, Ek,\so = %Ek,h- The increase in Ek.h at 
the boundaries of the convection zones are due to the 
horizontal deflection of the large sc ale flow and wave 
motions exc ited in stable layers (e.g. iMeakin fc Arnett 
20061 l2007al) and should be corrected for when identi- 
fying Ek,iso with the convective turbulence. 

In Fig. 3 (right) we over plot two approximations to 



Ek, iso'- one based on a uniform distribution of dissi- 
pation and one based on a dissipation that decreases 
linearly with enclosed mass. The relationship ejf = 
(2i?if i i so ) 3 / 2 /l c i with Id = H p is used. The absolute scale 
of the dissipation and kinetic energy profiles are pro- 
vided by the constraint that the global dissipation rate 
must balance the global rate of buoyancy driving. The 
amplitude of the kinetic energy that satisfies this global 
balance is found by varying it until the boundary condi- 
tions on Lk are satisfied (i.e., Lk = at the boundaries 
of the convection zone). Ek^so ~ \{Fcl P) 2 ^ 3 provides 
a good first approximation. 

The kinetic energy fluxes found using this procedure 
are compared to the simulation data in Fig. 4 for the 
two assumed dissipation profiles. 



in a future publication. For now we shall suffice to 
say that the dissipation can be derived directly from 
the stellar model by adopting certain constraints on 
the topology of the convective flow. In particular, a 
two component flow model consisting of a background 
isotropic turbulent state and a large scale, plume-like 
flow is a promising approach. 

The data presented in Figs. 3 and 4 illustrate the 
shortcomings of the commonly used closure relation re- 
ferred to as the down gradient approximation Most il- 
lustrative is the fact that while the Ek distributions are 
nearly identical in models hi and cl the Lk profiles are 
roughly mirror images. The locally defined down gra- 
dient approximation flux fails because the properties of 
the turbulent transport are strongly shaped by global 
constraints, a feature that is captured by the analysis 
presented in fjl] 

Another consequence of the distribution of kinetic 
energy within the convection zone is the rate of bound- 
ary layer mixing (see Fig. 2), which can significantly 
modify th e stellar structure on evo lutionary timescales 
(see §7 in lMeakin fc Arnettl[2007r} ). 
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5 Discussion 



We have provided a basic overview of the connection 
between turbulent dissipation and the kinetic energy 
flux in efficient (high Peclet number) convection. The 
only assumption made in our analysis involved the ra- 
dial profile of the dissipation ek which we will discuss 



1 The down gradient approximation is a closure relationship which 
relates the kinetic energy flux to the gradient of the turbulent 
kinetic energy such that Lk oc —X?Ek. 
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